function [ber_freq_div, fig_handle] = frequency_diversity(EbN0_dB, num_bits, num_subcarriers)
% 频率分集性能分析
% 输入:
%   EbN0_dB - 比特信噪比范围 (dB)
%   num_bits - 仿真比特数
%   num_subcarriers - 子载波数量数组
% 输出:
%   ber_freq_div - 频率分集误比特率矩阵
%   fig_handle - 图形句柄

    EbN0_linear = 10.^(EbN0_dB/10);
    ber_freq_div = zeros(length(EbN0_dB), length(num_subcarriers));
    
    % 添加路径
    addpath('../Common');
    colors = color_definitions();
    
    fprintf('=== 频率分集性能分析 ===\n');
    
    % 系统参数
    bandwidth = 1e6; % 1 MHz带宽
    coherence_bandwidth = 200e3; % 200 kHz相干带宽
    frequency_correlation = exp(-bandwidth / coherence_bandwidth); % 频率相关性
    
    for sub_idx = 1:length(num_subcarriers)
        N = num_subcarriers(sub_idx);
        fprintf('仿真 %d 子载波频率分集...\n', N);
        
        for eb_idx = 1:length(EbN0_dB)
            fprintf('  Eb/N0 = %d dB...\n', EbN0_dB(eb_idx));
            
            % 生成随机比特
            tx_bits = randi([0, 1], num_bits, 1);
            
            % BPSK调制
            tx_symbols = 2 * tx_bits - 1;
            
            % 频率选择性衰落信道
            % 生成相关频率响应
            freq_response = zeros(num_bits, N);
            for bit_idx = 1:num_bits
                % 生成频率响应 (相关瑞利衰落)
                base_response = sqrt(0.5) * (randn(N, 1) + 1i * randn(N, 1));
                
                % 应用频率相关性
                correlation_matrix = frequency_correlation .^ (abs((1:N)' - (1:N)));
                L = chol(correlation_matrix, 'lower');
                correlated_response = L * base_response;
                
                freq_response(bit_idx, :) = correlated_response';
            end
            
            % 添加AWGN噪声
            noise_power = 1 / (2 * N * EbN0_linear(eb_idx)); % 总功率分配
            noise = sqrt(noise_power) * (randn(num_bits, N) + 1i * randn(num_bits, N));
            
            % 接收信号
            rx_symbols = tx_symbols * ones(1, N) .* freq_response + noise;
            
            % 合并技术
            % 1. 选择合并 (选择最佳子载波)
            [~, best_subcarrier] = max(abs(freq_response), [], 2);
            rx_sc = zeros(num_bits, 1);
            for bit_idx = 1:num_bits
                rx_sc(bit_idx) = rx_symbols(bit_idx, best_subcarrier(bit_idx));
            end
            rx_bits_sc = (real(rx_sc) > 0);
            ber_sc = sum(tx_bits ~= rx_bits_sc) / num_bits;
            
            % 2. 等增益合并 (EGC)
            rx_egc = sum(rx_symbols, 2);
            rx_bits_egc = (real(rx_egc) > 0);
            ber_egc = sum(tx_bits ~= rx_bits_egc) / num_bits;
            
            % 3. 最大比合并 (MRC)
            weights = conj(freq_response) ./ (abs(freq_response).^2 + noise_power);
            rx_mrc = sum(weights .* rx_symbols, 2);
            rx_bits_mrc = (real(rx_mrc) > 0);
            ber_mrc = sum(tx_bits ~= rx_bits_mrc) / num_bits;
            
            % 保存MRC结果
            ber_freq_div(eb_idx, sub_idx) = ber_mrc;
        end
    end
    
    % 绘制结果
    fig_handle = figure('Name', '频率分集性能', 'Position', [400, 400, 1200, 800]);
    
    % 频率分集性能
    subplot(2,3,1);
    for sub_idx = 1:length(num_subcarriers)
        semilogy(EbN0_dB, ber_freq_div(:, sub_idx), ...
                 'Color', colors(sub_idx,:), 'LineWidth', 2);
        hold on;
        
        % 理论曲线
        N = num_subcarriers(sub_idx);
        theory_ber = nchoosek(2*N-1, N) * (0.5 * (1 - sqrt(EbN0_linear ./ (1 + EbN0_linear)))).^N;
        semilogy(EbN0_dB, theory_ber, 'Color', colors(sub_idx,:), 'LineStyle', '--', 'LineWidth', 1);
    end
    grid on;
    xlabel('Eb/N0 (dB)');
    ylabel('误比特率 (BER)');
    title('频率分集性能 (MRC合并)');
    legend([arrayfun(@(x) sprintf('N=%d', x), num_subcarriers, 'UniformOutput', false), ...
           arrayfun(@(x) sprintf('N=%d(理论)', x), num_subcarriers, 'UniformOutput', false)]);
    
    % 合并技术比较
    subplot(2,3,2);
    N = 8; % 固定8子载波
    ber_methods = zeros(length(EbN0_dB), 3);
    
    for eb_idx = 1:length(EbN0_dB)
        tx_bits = randi([0, 1], num_bits, 1);
        tx_symbols = 2 * tx_bits - 1;
        
        % 生成频率响应
        freq_response = zeros(num_bits, N);
        for bit_idx = 1:num_bits
            base_response = sqrt(0.5) * (randn(N, 1) + 1i * randn(N, 1));
            correlation_matrix = frequency_correlation .^ (abs((1:N)' - (1:N)));
            L = chol(correlation_matrix, 'lower');
            correlated_response = L * base_response;
            freq_response(bit_idx, :) = correlated_response';
        end
        
        noise_power = 1 / (2 * N * EbN0_linear(eb_idx));
        noise = sqrt(noise_power) * (randn(num_bits, N) + 1i * randn(num_bits, N));
        rx_symbols = tx_symbols * ones(1, N) .* freq_response + noise;
        
        % 选择合并
        [~, best_subcarrier] = max(abs(freq_response), [], 2);
        rx_sc = zeros(num_bits, 1);
        for bit_idx = 1:num_bits
            rx_sc(bit_idx) = rx_symbols(bit_idx, best_subcarrier(bit_idx));
        end
        ber_methods(eb_idx, 1) = sum(tx_bits ~= (real(rx_sc) > 0)) / num_bits;
        
        % 等增益合并
        rx_egc = sum(rx_symbols, 2);
        ber_methods(eb_idx, 2) = sum(tx_bits ~= (real(rx_egc) > 0)) / num_bits;
        
        % 最大比合并
        weights = conj(freq_response) ./ (abs(freq_response).^2 + noise_power);
        rx_mrc = sum(weights .* rx_symbols, 2);
        ber_methods(eb_idx, 3) = sum(tx_bits ~= (real(rx_mrc) > 0)) / num_bits;
    end
    
    methods = {'选择合并', '等增益合并', '最大比合并'};
    for method_idx = 1:3
        semilogy(EbN0_dB, ber_methods(:, method_idx), ...
                 'Color', colors(method_idx,:), 'LineWidth', 2);
        hold on;
    end
    grid on;
    xlabel('Eb/N0 (dB)');
    ylabel('误比特率 (BER)');
    title('合并技术比较 (N=8)');
    legend(methods);
    
    % 频率选择性衰落特性
    subplot(2,3,3);
    % 生成频率响应示例
    N_example = 64;
    freq_response_example = zeros(100, N_example);
    for i = 1:100
        base_response = sqrt(0.5) * (randn(N_example, 1) + 1i * randn(N_example, 1));
        correlation_matrix = frequency_correlation .^ (abs((1:N_example)' - (1:N_example)));
        L = chol(correlation_matrix, 'lower');
        freq_response_example(i, :) = abs((L * base_response)');
    end
    
    % 绘制频率响应
    freq_axis = (1:N_example) / N_example;
    mean_response = mean(freq_response_example, 1);
    std_response = std(freq_response_example, 1);
    
    plot(freq_axis, mean_response, 'b-', 'LineWidth', 2);
    hold on;
    plot(freq_axis, mean_response + std_response, 'r--', 'LineWidth', 1);
    plot(freq_axis, mean_response - std_response, 'r--', 'LineWidth', 1);
    grid on;
    xlabel('归一化频率');
    ylabel('幅度响应');
    title('频率选择性衰落特性');
    legend('均值', '±1标准差');
    
    % 相干带宽影响
    subplot(2,3,4);
    coherence_bandwidths = [50e3, 100e3, 200e3, 500e3]; % Hz
    N = 16; % 固定子载波数
    ber_coherence = zeros(length(EbN0_dB), length(coherence_bandwidths));
    
    for cb_idx = 1:length(coherence_bandwidths)
        Bc = coherence_bandwidths(cb_idx);
        freq_corr = exp(-bandwidth / Bc);
        
        for eb_idx = 1:length(EbN0_dB)
            tx_bits = randi([0, 1], num_bits, 1);
            tx_symbols = 2 * tx_bits - 1;
            
            freq_response = zeros(num_bits, N);
            for bit_idx = 1:num_bits
                base_response = sqrt(0.5) * (randn(N, 1) + 1i * randn(N, 1));
                correlation_matrix = freq_corr .^ (abs((1:N)' - (1:N)));
                L = chol(correlation_matrix, 'lower');
                correlated_response = L * base_response;
                freq_response(bit_idx, :) = correlated_response';
            end
            
            noise_power = 1 / (2 * N * EbN0_linear(eb_idx));
            noise = sqrt(noise_power) * (randn(num_bits, N) + 1i * randn(num_bits, N));
            rx_symbols = tx_symbols * ones(1, N) .* freq_response + noise;
            
            % MRC合并
            weights = conj(freq_response) ./ (abs(freq_response).^2 + noise_power);
            rx_mrc = sum(weights .* rx_symbols, 2);
            rx_bits = (real(rx_mrc) > 0);
            
            ber_coherence(eb_idx, cb_idx) = sum(tx_bits ~= rx_bits) / num_bits;
        end
    end
    
    for cb_idx = 1:length(coherence_bandwidths)
        semilogy(EbN0_dB, ber_coherence(:, cb_idx), ...
                 'Color', colors(cb_idx,:), 'LineWidth', 2);
        hold on;
    end
    grid on;
    xlabel('Eb/N0 (dB)');
    ylabel('误比特率 (BER)');
    title('相干带宽影响 (N=16)');
    legend(arrayfun(@(x) sprintf('Bc=%.0fkHz', x/1e3), coherence_bandwidths, 'UniformOutput', false));
    
    % OFDM系统仿真
    subplot(2,3,5);
    % OFDM参数
    N_fft = 64; % FFT大小
    N_cp = 16; % 循环前缀长度
    num_symbols = 100; % OFDM符号数
    
    % 频率选择性信道
    channel_taps = 8;
    channel_delays = (0:channel_taps-1) * (N_cp / channel_taps);
    channel_gains = sqrt(0.5) * (randn(1, channel_taps) + 1i * randn(1, channel_taps));
    
    ber_ofdm = zeros(1, length(EbN0_dB));
    
    for eb_idx = 1:length(EbN0_dB)
        total_errors = 0;
        total_bits = 0;
        
        for sym_idx = 1:num_symbols
            % 生成数据
            data_bits = randi([0, 1], N_fft, 1);
            data_symbols = 2 * data_bits - 1;
            
            % IFFT
            time_signal = ifft(data_symbols, N_fft);
            
            % 添加循环前缀
            tx_signal = [time_signal(N_fft-N_cp+1:end); time_signal];
            
            % 通过频率选择性信道
            rx_signal = zeros(N_fft + N_cp, 1);
            for tap = 1:channel_taps
                if channel_delays(tap) < length(rx_signal)
                    start_idx = channel_delays(tap) + 1;
                    end_idx = min(start_idx + length(tx_signal) - 1, length(rx_signal));
                    sig_start = 1;
                    sig_end = end_idx - start_idx + 1;
                    rx_signal(start_idx:end_idx) = rx_signal(start_idx:end_idx) + ...
                        channel_gains(tap) * tx_signal(sig_start:sig_end);
                end
            end
            
            % 添加噪声
            noise_power = 1 / (2 * EbN0_linear(eb_idx));
            noise = sqrt(noise_power) * (randn(N_fft + N_cp, 1) + 1i * randn(N_fft + N_cp, 1));
            rx_signal = rx_signal + noise;
            
            % 移除循环前缀
            rx_no_cp = rx_signal(N_cp+1:end);
            
            % FFT
            rx_symbols = fft(rx_no_cp, N_fft);
            
            % 均衡 (假设已知信道)
            channel_freq = fft([channel_gains, zeros(1, N_fft-channel_taps)], N_fft);
            rx_equalized = rx_symbols ./ channel_freq;
            
            % 判决
            rx_bits = (real(rx_equalized) > 0);
            
            % 计算误码
            total_errors = total_errors + sum(data_bits ~= rx_bits);
            total_bits = total_bits + N_fft;
        end
        
        % 确保total_errors和total_bits是标量
        total_errors = sum(total_errors(:));
        total_bits = sum(total_bits(:));
        ber_ofdm(eb_idx) = total_errors / total_bits;
    end
    
    semilogy(EbN0_dB, ber_ofdm, 'o-', 'LineWidth', 2, 'Color', colors(1,:));
    hold on;
    
    % 理论OFDM性能 (AWGN)
    theory_awgn = 0.5 * erfc(sqrt(EbN0_linear));
    semilogy(EbN0_dB, theory_awgn, 'k--', 'LineWidth', 2);
    
    grid on;
    xlabel('Eb/N0 (dB)');
    ylabel('误比特率 (BER)');
    title('OFDM系统性能');
    legend('OFDM仿真', 'AWGN理论');
    
    % 子载波分配策略
    subplot(2,3,6);
    allocation_strategies = {'随机分配', '最佳选择', '注水分配'};
    N = 16;
    ber_allocation = zeros(length(EbN0_dB), length(allocation_strategies));
    
    for alloc_idx = 1:length(allocation_strategies)
        for eb_idx = 1:length(EbN0_dB)
            tx_bits = randi([0, 1], num_bits, 1);
            tx_symbols = 2 * tx_bits - 1;
            
            % 生成频率响应
            freq_response = zeros(num_bits, N);
            for bit_idx = 1:num_bits
                base_response = sqrt(0.5) * (randn(N, 1) + 1i * randn(N, 1));
                correlation_matrix = frequency_correlation .^ (abs((1:N)' - (1:N)));
                L = chol(correlation_matrix, 'lower');
                correlated_response = L * base_response;
                freq_response(bit_idx, :) = correlated_response';
            end
            
            % 子载波分配
            if alloc_idx == 1 % 随机分配
                selected_subcarriers = randi(N, num_bits, 1);
            elseif alloc_idx == 2 % 最佳选择
                [~, selected_subcarriers] = max(abs(freq_response), [], 2);
            else % 注水分配
                % 简化的注水分配
                avg_gains = mean(abs(freq_response).^2, 1);
                [~, sorted_idx] = sort(avg_gains, 'descend');
                selected_subcarriers = mod((1:num_bits)', length(sorted_idx)) + 1;
                selected_subcarriers = sorted_idx(selected_subcarriers);
            end
            
            % 分配子载波
            allocated_response = zeros(num_bits, 1);
            for bit_idx = 1:num_bits
                allocated_response(bit_idx) = freq_response(bit_idx, selected_subcarriers(bit_idx));
            end
            
            % 添加噪声
            noise_power = 1 / (2 * EbN0_linear(eb_idx));
            noise = sqrt(noise_power) * (randn(num_bits, 1) + 1i * randn(num_bits, 1));
            rx_symbols = tx_symbols .* allocated_response + noise;
            
            % 判决
            rx_bits = (real(rx_symbols) > 0);
            ber_allocation(eb_idx, alloc_idx) = sum(tx_bits ~= rx_bits) / num_bits;
        end
    end
    
    for alloc_idx = 1:length(allocation_strategies)
        semilogy(EbN0_dB, ber_allocation(:, alloc_idx), ...
                 'Color', colors(alloc_idx,:), 'LineWidth', 2);
        hold on;
    end
    grid on;
    xlabel('Eb/N0 (dB)');
    ylabel('误比特率 (BER)');
    title('子载波分配策略 (N=16)');
    legend(allocation_strategies);
    
    % 保存图片
    saveas(fig_handle, 'results/frequency_diversity_analysis.png');
end